﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MathLib {
    public class GenerationUtils {


        //Вектор начального состояния популяции
        public static Matrix GenerateInit(double popValue, int groupCount) {
            double last = Math.Sqrt(popValue);
            Gauss g = new Gauss();
            Matrix v = new Matrix(groupCount, 1);
            int first = 0;
            double value;
            double h = (last - first) / groupCount;
            double step = 0;
            for (int i = 0; i < groupCount; i++) {
                value = g.NextGaussian(0, last / 100);
                v[i, 0] = Math.Round(-(step + value) * (step + value) + (last + value) * (last + value));
                step += h;
            }
            return v;
        }


        public static Matrix GenerateAvgPopValue(int groupCount, Matrix[] populations) {
            Matrix v = new Matrix(groupCount, 1); ;

            for (int j = 0; j < groupCount; j++) {
                double s = 0;
                for (int i = 0; i < populations.Length; i++) {
                    s = s + populations[i][j, 0];
                }
                v[j, 0] = s / populations.Length;
            }
            return v;
        }

        //Вектор шума
        public static Matrix[] GenerateNoise(int groupCount, Matrix[] populations, double sigma, bool fullind) {

            Matrix[] noise = new Matrix[populations.Length];
            int n = populations.Length;

            if (sigma == 0) {
                for (int j = 0; j < n; j++) {
                    noise[j] = new Matrix(groupCount, 1);
                }
                return noise;
            }

            Gauss g = new Gauss();


            if (fullind) {
                for (int j = 0; j < n; j++) {
                    noise[j] = new Matrix(groupCount, 1);
                    for (int i = 0; i < groupCount; i++) {
                        double stdDev = populations[j][i, 0] * sigma;
                        double value = g.NextGaussian(0, stdDev);
                        noise[j][i, 0] = value;
                    }
                }
                return noise;
            } else {
                
                for (int j = 0; j < n; j++) {
                    noise[j] = new Matrix(1, 1);
                    double s = 0;
                    for (int i = 0; i < populations[0].NoRows; i++) {
                        s += populations[j][i, 0];
                    }
                    double stdDev = s * sigma;
                    double value = g.NextGaussian(0, stdDev);
                    noise[j][0, 0] = value;
                }    
            }
            return noise;
        }


        public static Matrix[] GenerateMatrixNoise(int groupCount, Matrix[] populationMatrixes, double utf, double utm) {

            Matrix[] noise = new Matrix[populationMatrixes.Length];

            if (utm == 0 && utf == 0) {
                for (int j = 0; j < populationMatrixes.Length; j++) {
                    noise[j] = new Matrix(groupCount, groupCount);
                }
                return noise;
            }

            Gauss g = new Gauss();
            double stdDev;

            for (int k = 0; k < populationMatrixes.Length; k++) {
                noise[k] = new Matrix(groupCount, groupCount);

                for (int i = 0; i < groupCount; i++) {
                    stdDev = populationMatrixes[k][0, i] * utf;
                    double value = g.NextGaussian(0, stdDev);
                    noise[k][0, i] = value;
                }

                for (int i = 1; i < groupCount; i++)
                    for (int j = 0; j < groupCount; j++) 
                     if (i - 1 == j)
                     {
                        stdDev = populationMatrixes[k][i, j] * utm;
                        double value = g.NextGaussian(0, stdDev);
                        noise[k][i, j] = value;
                    }

            }

            return noise;
        }

        

        //Вектор рождаемости
        private static Matrix GenerateFertilityRate(int groupCount, double maxFertility) {
            Gauss g = new Gauss();
            Matrix v = new Matrix(groupCount, 1);
            int first = -groupCount;
            int last = groupCount;
            double h = (last - first) / groupCount;
            double value;
            double x = first;
            for (int i = 0; i < groupCount; i++) {
                value = g.NextGaussian(0, last / 100);
                v[i, 0] = Math.Round(maxFertility * Math.Exp(-(x * x) / 10) + Math.Abs(value), 4);
                x = x + h;
            }
            return v;
        }

        //Вектор смертности
        private static Matrix GenerateMortalityRate(int groupCount) {
            Gauss g = new Gauss();
            Matrix v = new Matrix(groupCount - 1, 1);
            double value;
            for (int i = 0; i < groupCount - 1; i++) {
                value = g.NextGaussian(0, 0.01);
                v[i, 0] = Math.Round(1 - 0.113 * Math.Sqrt(i) - Math.Abs(value), 4);
            }
            return v;
        }

        //Матрица A
        public static PopMat GenerateA(int groupCount, double maxFertility) {
            PopMat popMat = new PopMat();
            Matrix f = GenerateFertilityRate(groupCount, maxFertility);
            Matrix m = GenerateMortalityRate(groupCount);
            MatrixA A = new MatrixA(groupCount, f, m);
            popMat.A = A;
            popMat.matrixMainValue = A.IsStable;
           return popMat;
        }


        //Матрица A
        public static PopMat GenerateAReal(int groupCount, Matrix f, Matrix m) {
            PopMat popMat = new PopMat();
            MatrixA A = new MatrixA(groupCount, f, m);
            popMat.A = A;
            popMat.matrixMainValue = A.IsStable;
            return popMat;
        }


        private static double CheckStability(Matrix f, Matrix m) {
            double c = 1;
            double r = 1;
            for (int i = 0; i < f.NoRows; i++) {
                c = 1;
                for (int j = 0; j < i - 1; j++)
                    c *= m[j, 0];
                r = f[i, 0] * c;

            }
            return r;
        }

        //Генерирует последовательность векторов состояний популяции
        public static GeneratedPopulation GeneratePopulation(int populationsCount, double maxFertility, Matrix p0) {

            Matrix[] PopulationMatrixes = new Matrix[populationsCount];
            Matrix[] populations = new Matrix[populationsCount];
            PopMat popmat = new PopMat();
            popmat = GenerateA(p0.NoRows, maxFertility);

            populations[0] = p0;
            PopulationMatrixes[0] = popmat.A;
            for (int i = 1; i < populationsCount; i++) {
                populations[i] = popmat.A * populations[i - 1];
                popmat = GenerateA(p0.NoRows, maxFertility);
                PopulationMatrixes[i] = popmat.A;
            }

            GeneratedPopulation pop = new GeneratedPopulation();
            pop.populationMatrixes = PopulationMatrixes;
            pop.populations = populations;
            pop.matrixMainValue = popmat.matrixMainValue;
            return pop;
        }



        //Генерирует последовательность векторов состояний реальной популяции
        public static GeneratedPopulation RealPopulation(int populationsCount, Matrix p0, Matrix f, Matrix m) {

            Matrix[] PopulationMatrixes = new Matrix[populationsCount];
            Matrix[] populations = new Matrix[populationsCount];
            PopMat popmat = new PopMat();
            popmat = GenerateAReal(p0.NoRows, f,m);

            populations[0] = p0;
            PopulationMatrixes[0] = popmat.A;
            for (int i = 1; i < populationsCount; i++) {
                populations[i] = popmat.A * populations[i - 1];
                popmat = GenerateAReal(p0.NoRows, f, m);
                PopulationMatrixes[i] = popmat.A;
            }

            GeneratedPopulation pop = new GeneratedPopulation();
            pop.populationMatrixes = PopulationMatrixes;
            pop.populations = populations;
            pop.matrixMainValue = popmat.matrixMainValue;
            return pop;
        }



        //Матрица наблюдения
        public static Matrix GenerateC(int groupCount, bool fullind) {

            if (fullind) { //Единичная матрица наблюдения
                Matrix C = new Matrix(groupCount, groupCount);
                for (int i = 0; i < groupCount; i++)
                    for (int j = 0; j < groupCount; j++)
                        C[i, j] = i == j ? 1.0 : 0.0;
                return C;
            } else 
            {
                Matrix C = new Matrix(groupCount, 1); //Вектор-строка наблюдения
                for (int i = 0; i < groupCount; i++)
                        C[i,0] = 1;
                return C;
            }
            
        }


        // Матрица ковариации ошибки оценки начального состояния
        public static Matrix GenerateP(int groupCount, double P0) {
            Gauss g = new Gauss();
            Matrix m = new Matrix(groupCount, groupCount);
            for (int i = 0; i < m.NoRows; i++)
                for (int j = 0; j < m.NoRows; j++)
                    if (i == j) m[i, j] = Math.Abs(g.NextGaussian(P0, 0));
                    else m[i, j] = 0;
            return m;
        }



        //Матрица ковариации шумов наблюдения
        public static Matrix GenerateVnu(int groupCount, double vnu, Matrix avgPopulationValue, bool fullind) {

            if (fullind) {

                Matrix m = new Matrix(groupCount, groupCount);

                for (int i = 0; i < m.NoRows; i++)
                    for (int j = 0; j < m.NoRows; j++)
                        if (i == j) m[i, j] = vnu * vnu * avgPopulationValue[j, 0] * avgPopulationValue[j, 0];
                        else m[i, j] = 0;
                return m;
            } else {
                Matrix m = new Matrix(1, 1);
                double s = 0;
                for (int j = 0; j < avgPopulationValue.NoRows; j++)
                    s += avgPopulationValue[j, 0];
                m[0, 0] = vnu * vnu * s * s;
                return m;
            }
        }

        //Матрица ковариации ошибок рождаемости и передвижки
        public static Matrix GenerateUt(int groupCount, double utf, double utm, Matrix avgPopulationValue, Matrix populationMatrixes) {

            Matrix m = new Matrix(groupCount, groupCount);
            double value = 0 ;

            for (int i = 0; i < m.NoRows; i++)
                for (int j = 0; j < m.NoRows; j++) 
                {
                    if (i == 0) value += utf * utf * avgPopulationValue[j, 0] * avgPopulationValue[j, 0] * populationMatrixes[0, j] * populationMatrixes[0, j];
                    if (i == j && i != 0) m[i,j] =  utm * utm * avgPopulationValue[j, 0] * avgPopulationValue[j,0] * populationMatrixes[i-1, j] * populationMatrixes[i-1, j];

                    else m[i, j] = 0;
                }
            m[0,0] = value;
            return m;
        }

        public static Matrix GenerateDisp(int groupCount, Matrix AvgPopValue, double vnu) {
            Matrix d = new Matrix(groupCount, 1);
            for (int i = 0; i < d.NoRows; i++)
                d[i,0]= AvgPopValue[i,0]*AvgPopValue[i,0]*vnu*vnu;
                    return d;
        }

        public static Matrix MortalityEstimation(Matrix x_prev, Matrix x_next, Matrix A, int groupCount) {
            Matrix newA = new Matrix(groupCount, groupCount);

            for (int i = 0; i < groupCount; i++)
                for (int j = 0; j < groupCount; j++)
                    newA[i, j] = A[i, j];

            Matrix newm = new Matrix(groupCount - 1, 1);

            for (int i = 1; i < groupCount; i++) 
                newm[i - 1, 0] = x_next[i, 0] / x_prev[i - 1, 0];
       

            for (int i = 1; i < groupCount; i++)
                newA[i, i - 1] = newm[i - 1, 0];

            return newA;
        }



    }



    public struct PopMat {
        public double matrixMainValue;
        public Matrix A;
    }

    public struct GeneratedPopulation {
        public Matrix[] populationMatrixes;
        public Matrix[] populations;
        public double matrixMainValue;
    }

}
